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of LAN.DSAT images is demonstrated. 


This report describes research done at the Artificial Intelligence Laboratory 
of the Massachusetts Institute of Technology. Support for the laboratory's 
artificial intelligence research is provided in part by the Advanced Research 
Projects Agency of the Department of Defense under Office of Naval Research 
contract N00014-75-C-0643. 


@ Massachusetts Institute of 
Technology 1978 








- 1 - 


1. Destripinq of images obtained using multiple sensors. 

An image sensing device using a single photoelectric sensor which is 
mechanically scanned across the scene produces outstanding digitized images 
since sensitivity, resolution and transfer functions are the same for all 
points in the image. Unfortunately, such a device is limited in speed by the 
mechanical movement. More importantly, it is limited in speed by the fact 
that an accurate measurement of scene radiance requires the collection of an 
adequate number of photons. This explains the preponderance of linear arrays 
of sensors and area sensors such as vidicons which are otherwise deficient 
because of geometric distortions, non-uniform response, non-uniform resolution, 
and so on. 

A compromise can be struck, where a small set of sensors is mechanically 
scanned to collect the image. In the system used aboard LANDSAT, for example, 
each spectral band is scanned using six sensors at the same time. Thus, six 
lines of the image are produced during a single sweep of the mirror. On the 
next sweep the satellite has advanced in its orbit by an amount which allows 
the same set of sensors to pick up the next six lines of the image [1]. 

Unfortunately, the sensors do not have identical transfer functions. As 
a result, images produced in this fashion show undesirable, regular "striping". 
This effect can be removed if the transfer functions are accurately known, 
since one could then compute scene radiance from the sensor output using 
the inverses of these transfer functions. The sensors used in older equipment 
in particular have time-varying behavior. Photomultipliers, for example, show 
a drift in both gain and offset (dark current) due to small changes in the 
material of the dynodes used in the electron multiplier stages and temperature 
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variations. 

If a reference object containing all scene radiances of interest were 
in the scene, one could recalibrate the sensors continuously. This is diffi¬ 
cult to arrange. An alternative is the scanning of a gray wedge placed over 
a light source at the end of every scan line. This, in fact, is what is done 
aboard LANDSAT. The results are used to estimate the gains and offsets of the 
sensors. The digital data produced from the raw satellite signals is corrected 
using this information [1]. 

Unfortunately, one finds that the striping effect is not removed in this 
fashion; the reasons for this are not entirely clear. One cause appears to 
be the use of the calibration data as a means of adjusting gain and offset 
so that each sensor is related to its preflight condition. Slight changes in 
the light source, the gray wedge and the geometry of imaging introduce drifts 
which are not compensated for. Another reason is related to the fact that 
photomultipliers are somewhat nonlinear and have a response which depends on 
their exposure history. Modern devices using solid state photodiodes do not 
suffer from these problems. 

The methods explored here for destriping images are based on the assump¬ 
tion that each sensor is exposed to scene radiances with approximately the 
same probability distribution. The sensor values can then be modified so that 
each one is related in the same way to the actual scene radiance. The in¬ 
formation required to perform this modification is extracted from statistics 
of the observed sensor outputs. 






2. A simple method for linear transducers. 


If the image sensors are linear and time invariant, a simple method can 
be used to reduce striping. The sensor output, x', can be written as a func¬ 
tion of the scene radiance, x, as follows: 

x 1 = f(x) = a + b • x 

Each sensor has its own, fixed values of offset, a, and gain, b. If these 
are known, the scene radiance can be calculated using the inverse of the 
transfer function. 


x = g(x') = (x‘ - a)/b 

If this is done for each sensor in turn, striping effects will be removed. 

The required constants for each sensor can be determined if a calibration 
object containing two or more known scene radiance values is available in the 
scanned scene. If such a calibration object is not available one can estimate 
the (relative) values of gain and offset using simple statistics of observed 
sensor values. Each sensor sees a subimage consisting of every nth line 
(when n sensors are used). The complete image is formed by interlacing these 
subimages. It seems reasonable to suppose that, for a large enough image, 
each subimage has approximately the same probability distribution of scene 
radiance values. One would not expect a particular subimage to contain many 
more values in a particular range of scene radiances than another subimage. 

If this assumption is correct, then the gain and offset constants can be 






estimated from the mean and standard deviation of the measured sensor output 
values. If the mean of the scene radiance is y and the standard deviation is 
a, then the mean of the sensor output will be y ' = a + b y and the standard 
deviation of the sensor output a' = b a. Then, 

b = a 1 /a 

and 

a = (y 1 a - y o' ) /a 

Clearly, it is not reasonable to assume that one can find the absolute values 
of the mean and standard deviation of the actual scene radiance. Fortunately, 
for destriping purposes only relative values are important. That is, one can 
use the mean and standard deviation of the sensor outputs for the whole image 
in place of the mean and standard deviation of the scene radiance. Naturally 
now the results will not be scene radiance values. The striping however will 
be removed since each subimage now has the same mean and standard deviation, 
and, if the assumption introduced earlier applies, the same linear relation¬ 
ship to scene radiance. 

Note that one can relax the assumption about the relationship of the 
subimages. Here it is not necessary that they have the same probability dis¬ 
tribution of scene radiance, only that their means and standard deviations be 
the same. This simple method has been applied by some users of LANDSAT data 





3. Shortcomings of the simple method. 

We have found this method to be only partially successful in destriping 
LANDSAT images. One reason for this may be that out of a range of 128 possible 
sensor outputs a range of only around 30 values correspond to normal scene 
radiance values. Low values are not found in short wavelength bands because 
of light scatter in the air. Conversely, large values correspond to cloud, 
snow and ice, and scene radiance values of such areas often exceed the highest 
available sensor output values and so result in clipping . Clipping of sensor 
values corresponding to low scene radiances also occurs at times in the long 
wavelength bands due to negative sensor offsets. Both of these nonlinear ef¬ 
fects will introduce skew into the calculation of means and standard deviations. 

One may alleviate this problem by removing sensor values outside a certain 
range from consideration. While slightly better results are obtained in this 
fashion, it is clear that the arbitrarily selected thresholds needed introduce 
biases of their own. This later effect can be dealt with by eliminating the 
same fraction of sensor values from the low end of the output of each sensor. 
Similarly, a fixed fraction of sensor values is removed from the high end. 

Even with this refinement, results are not entirely satisfactory. Super¬ 
ficially, it appears that different gains and offsets are appropriate for 
different scene radiance ranges. That is, the sensor transfer curves are some¬ 
what nonlinear. We thus devised a method which deals with this problem directly. 
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4. Preliminary considerations. 

Consider a random variable X with probability density function p(x). 
The function p(x) is non-negative and satisfies 



p(x) dx = 1 



The probability density function p(x) can be estimated from a large number N 
of observations of the random variable X. If n of these measurements fall in 
the interval [x, x + <5x], then n/N tends to p(x) sx as N becomes very large 
and 6x small (in a fashion which allows N Sx, and thus n, to become large 


also). 


The cumulative probability density function P(x) is defined as 


r x 

P(x) = / p(t) dt 

* — oo 


This function is monotonically non-decreasing since p(x) is non-negative. P(x) 
represents the probability that the random variable X takes on a value less 
than or equal to x. 

Now consider observing the random variable X by means of a transducer 
with transfer function f(x). Its output can be thought of as a new random 
variable X', say, with a probability density function p'(x*). This function 
is related to the probability density function p(x) of the original random 
variable X, in a fashion which depends on the transfer function x‘ = f(x). 

It is easiest to develop this relationship in terms of the cumulative distri- 
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bution functions P(x) and P'(x') where 

i / 


P 



p 1 (t) dt 


If x 1 lies in a range R' when x lies in the range R, then clearly. 


j p'(x') dx' = J p(x) dx 
R* R 




Now assume that f(x) is monotonically non-decreasing. Then the range 
x S x Q is mapped into the range x' - f(x ). 

P'[f(x)] = P(x) 


As a result one can determine the transfer function f(x) if the cumulative 
probability density functions P(x) and P'(x) are known and if the latter has 
an inverse. Then, 

f(x) = (P 1 )" 1 P(x) 


If P 1 is monotonically increasing, the required inverse will exist. Difficulties 
will be encountered only when P'(x') is constant over a certain range. That is, 
if P 1 (x') = c [and hence p'(x 1 ) = 0] for x 1 e [x{, X2]. Then, if P(x) = c, one 
can say only that f(x) e [x[, X 2 ]. 

There are two possible causes of this problem. First it may be that f(x) 



actually has a discontinuity. In this case, one correctly finds a jump from 
x i to X 2 in the solution. The other possibility is more serious. If p'(x') = 
0 because p(x) = 0 [where x 1 = f(x) as before], then the transfer function 
f(x) cannot be found in the specified range because, in essence, no inputs are 

available to test it in this range. The information to recover f(x) there 
is thus not available. 

Note, however, that if the inputs to the transducer are in fact character 
ized by the given probability density function, then our lack of knowledge 
of the transfer function in the specified range is of no consequence since 
no inputs fall in this range anyway. 

To calculate scene radiance from sensor values, we actually need the in¬ 
verse g(x') of the transfer function. This can be found just as easily. If, 

P'(x') = P [g(x')] 


Then 


g(x') = P" 1 p'(x') 

The same considerations regarding the existence of the inverse P" 1 apply here 
as those discussed regarding the existence of the inverse (P') _1 . All these 

special case problems are avoided if the cumulative probability distribution 
functions are monotonically increasing. 

The method shown here for finding the transfer function of a transducer 
(or its inverse) is based on the same analysis as that used to design a genera¬ 
tor of pseudo-random numbers with a desired probability distribution function 





P'(x) when a generator is available which produces pseudo-random numbers with 
known probability distribution function p(x) [4]. 

A graphical illustration of the relationships discussed may be found in 
Figure 1. The dotted line suggests how one may determine the transfer func¬ 
tion value, f(x), given a scene radiance value x. Conversely, the same dotted 
line may be followed in the reverse direction to find the value, g(x'), of the 
inverse function from a given value of the transducer output x 1 . 
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5. Estimation from a finite number of samples. 

To apply this method to determine the transfer function of a real trans¬ 
ducer, the cumulative probability density functions must be determined from a 
model of the underlying process generating the random variables or estimated 
from frequencies of observed occurrence using a finite number of samples. In 
the latter case an uncertainty will be found in the estimation of the proba¬ 
bilities which will be inversely proportional to the square root of the number 
of samples falling in a particular interval. 

In particular, the sample deviation of the estimate of the cumulative 
;probability P(x) obtained from a total of N measurements is 

aP = ^P(l - P)/N 

The following approximate analysis may be helpful: Let of be the uncertainty 

in estimating the transfer function, x' = f(x), resulting from the uncertainty aP 
in estimating P(x), then, 

P'[f(x) + of] = P(x) + aP 

If P (x ) is reasonably well behaved then the left-hand side can be expanded, 

p 1 [f(x)] + p'(x') of ~ P(x) + aP 

using the fact that p' is the first derivative of P'. Clearly then 


p 1 (x') of ~ aP 
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The uncertainty in estimating the transfer function x 1 = f(x) is thus highest 
where p'(x') is small. In fact, we have already shown that f(x) cannot be 
uniquely determined when p'(x') = 0. 

One can similarly show that 

p(x) og = aP 1 

where ag is the uncertainty in estimating the inverse transfer function, 
x = g(x‘), due to uncertainty a P in estimating P'(x'). Not surprisingly, un¬ 
certainty in estimating the inverse transfer function is highest where p(x) 
is small. 

The probability density functions p(x) and p'(x') are related by, 

P(x) = P'(x') and P'Oc') = p(x) 

when the indicated derivatives exist. 
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6- Transducer with discrete output values. 

Essentially the same method may be used if the transducer produces dis 
Crete outputs. Consider, for example, a case where the input range can be 
broken up into a number of intervals such that 

f(x) = i if x e [x., x f+1 ) 

The probability density function of the output of the transducer is then di 
Crete and, 

11m {* i+1 ‘ e 

P i = e + 0 I P(x) dx 

X i 

Clearly, pj > 0 and 

co 

I Pj = 1 

i = -°° 


The cumulative probability density function can be defined as follows 


1 



If f(x) is monotonically non-decreasing, then the same argument applied in 
the continuous case, leads again to 

p, -£f(*)]- p W 




If P 1 can be inverted, the transfer function can be found using 


f(x) = (P') _1 P(x) 

The only difference is that here f(x) is a function from a continuous range 
to a discrete domain. Naturally, when one finds the inverse of the transfer 
function, g(x'), using these methods, one has to accept the fact that the 
actual value of x cannot be recovered, only a range [x., x. ) 

1 1 + I 1 * 




-14- 


7. Application to satellite images, 

A particular sensor, in a system using n sensors, sees a subimage contain- 
ing every n line. To apply the methods developed here one has to assume that 
each sensor is exposed to scene radiances with similar probability distribu¬ 
tions. If the image is large enough, it is very unlikely that one sensor will 
see substantially fewer or more scene radiances in a particular range. A 
sensor's properties can then be estimated from the statistics of its outputs. 

Since the probability distribution function of the actual scene radiance 
is not available, only relative adjustments can be made. That is, the proba¬ 
bility distribution function of sensor outputs for the whole image is used as 
a reference. Consequently, application of the inverse functions so determined 
to each subimage will not produce scene radiances. It will, however, result 
in an image in which each image gray level is related to the scene radiance 
in the same way. Thus, striping will have been removed. 

A graphical illustration of the method is shown in Figure 2. Note that 
here both P(x) and P'(x') are discrete. This results in a small problem 
when the transfer function f(x) is to be found using data from a real image. 

For perfect dats, there always is a value x' for every value x, such that 
P'(x') = P(x). As indicated by the dotted line in the figure, this may not be 
the case when the data is obtained from a finite number of samples obtained 
by different sensors. The best one can do then is to find a value x' such that 


P'(x') s P(x) < P'(x' + 1) 


111 Similarly, in determining' the inverse transfer function g(x'), one can cfo no 
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better than finding a value x such that 
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P(x) < P'(x') < P(x + 1) 


What value should be used in the lookup table for destriping? One might argue 


that some values should be translated to x, others to(x +1). If this is done 
in the appropriate fashion, the histogram of gray levels in the destriped image 
will e qual the histogram of gray levels in the raw image. It is difficult 


however to decide which points should receive one value and which the other. 
If the selection is based on a random number generator, additional noise will 
be introduced. In any case a change smaller than one gray level is usually 
negligible considering the limitations of radiometric accuracy in the imaging 


system. The algorithm we employed, arbitrarily uses the smaller value of the 
two possible ones. Note that the histogram of the destriped image will not 
be e xactly equal to the histogram of the raw image when this is done. 
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8. Details of the algorithms. 

The first step is the determination of a cumulative histogram of sensor 
values for the whole image as a reference. Let there be H(x) occurrences of 
sensor outputs less than or equal to x out of a total of N values. Now for 
the subimage produced by sensor i, one calculates a similar cumulative histo¬ 
gram. Let H.(x') be the number of sensor outputs less than or equal to x 1 , 
produced by sensor i, out of a total of N. values. Here, 

n 

N = I N. 
i = 1 1 

where n is the number of sensors. 

A lookup table g(x') is now constructed by applying the inverse of the 
function H(x) to H^(x'). This lookup table is then used to modify all the 
sensor values produced by sensor i. The inverse can be calculated relatively 
easily since H(x) is a monotonically non-decreasing function. The lookup 
table value g(x') is the smallest number x such that 

N. H(x) * N H.j (x') < N. H(x + 1) 

This process is repeated for each sensor in turn, until all image values have 
been modified by the lookup table appropriate to the sensor with which they 
were measured. 






9. Results. 


The image used for the illustrations here is a part of LANDSAT-1 image 
E-l 078-09555 taken 1972/0ctober/9. The part used contains 364 lines each of 
430 pixels. The worst striping in this case occurred in band 6 (,7y to .8y) 
and band 7 (.8y to l.ly), so the discussion will concentrate on these two. 

The data transmitted from the spacecraft is pseudo-logarithmically compressed 
and converted to six bit numbers. On the ground the compression is removed 
using a lookup table. The result is a seven bit number. Naturally, about 
half of the possible seven bit numbers never occur in a given subimage. Fur¬ 
thermore, scene radiance corresponding to normal surfaces (Other than ice, 
snow or cloud) occupy only the lower 25% to 35% of the available range. As a 
result, the total range of scene radiances of interest corresponds to a very 

small set of distinct sensor output values. This contributes a little to diffi 
culties in completely removing striping effects. 

Comparison of Figure 3, the original band 6 image, with Figure 4, the pro¬ 
cessed version, shows that the heavy, regular striping is removed by the pro¬ 
cessing described here. It is instructive to inspect the lookup tables used 
for each sensor. These are shown as six subfigures in Figure 5. The short 
horizontal sections in the transfer function correspond to sensor values which 
never occur as a result of the decompression table lookup. The inverse trans¬ 
fer functions shown in Figure 5 appear to evidence some nonlinearities. This 

explains why the simple destriping technique described earlier does not do as 
wel 1. 

Comparison of Figure 6, the original band 7 image, with Figure 7, the 
processed version, shows that the regular striping is removed here as well. 






The very light striping remaining in shadow areas amounts to fluctuations of only 
one or two pixel values, but may be discernable because of the human observer's 
sensitivity to small differences of reflectance in dark areas. Inspection of the 
histograms of raw image data reveals uneven intervals in the analog-to-digital 
convertors used, which contribute to this effect (in fact some codes near major 
transitions never occur). The relative smoothness of the lookup tables used in 
destriping this band, shown in Figure 8, is in part due to the fact that no 
compression and decompression is performed on data for this band; six bit 
numbers linearly related to sensor output being transmitted directly. The 
inverse transfer functions shown in Figure 8 appear to be fairly linear, which 
is probably a result of the linearity of the silicon photodiodes used for this 
infrared band. One would expect then that the simple destriping method would 
be fairly successful for this band. 

We experienced some difficulties due to the normalization processing per¬ 
formed on the raw LANDSAT data. It may be useful to provide users of LANDSAT 
tapes optionally with the original data. Preliminary results indicate that 
slightly better destriping may be possible using the raw image sensor values. 

It is also unfortunate that areas of high reflectance produce scene radiances 
which saturate the imaging system. If this was not the case, image sensor 
outputs corresponding to areas covered by thick clouds could be used in cali¬ 
bration of absolute reflectances as well as in normalization for destriping 


purposes. 



10. Relation to Histogram Normalization Methods. 


A number of techniques aimed at the enhancement of images intended for 
human viewing are based on manipulation of the gray level histogram. Some, for 
example, transform the histogram into one considered more desirable [5, 6, 7, 

8, 9], either flat or "hyperbolic". The suggestion has been made that such 
techniques may allow for changes in sensor characteristics or scene illumination 
[6, 8, 10]. Several suggestions have been made regarding the difficulty in¬ 
dicated in Figure 2 relating to the mismatch of two cumulative histograms. 

One can, for example, increase the apparent fineness of quantization by taking 
into account the context of each picture cell [7]. If each picture cell has a 
maximum gray level value m and k neighbors, one may multiply its gray level 
value by mk and subtract (or add) the sum of the neighboring gray level values. 
In this way rank-ordering of picture cell values is preserved, while the number 
of possible gray level values is vastly increased, bringing the distributions 
closer to those found in the continuous case. 

Recent work on matching of images obtained by one sensor at different 
times is perhaps most closely related to our work here on destriping [10]. This 
sort of approach has however not yet found much of a following in the remote 
sensing community [11]. We believe our application of histogram "equalization" 
to subimages obtained using sensors of a multiple-sensor system is novel. 
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11. Conclusions. 

The following can be accomplished using the method presented here: 

1. The effect of differences in the transfer functions of the 
different sensors is removed. All gray levels are then 
related in the same way to the original scene radiance values. 

2. The overall tone-scale is preserved. That is, the histogram 
of gray levels of the destriped image is (approximately) 

the same as the histogram of the raw image. There is no loss 
in resolution, nor is the noise level increased. 

3. The computation requires only two serial passes over the whole 

image: one to collect the relevant histograms, another to 

apply the inverse transducer function represented as a simple 
lookup table. 
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13. Figures. 


Figure 1. 


Figure 2. 

Figure 3. 

Figure 4. 
Figure 5. 


Figure 6. 

Figure 7. 
Figure 8. 


The transfer function, x' = f(x), of a transducer can be found if 
the cumulative probability distribution functions of its input, x, 
and its output, x', are known. 

When both cumulative probability distribution functions are discrete 

there may not be a value of x' so that P(x') = P(x). Some rule must 
be adopted to deal with this. 

Original, unrectified image of band 6 (,7u to .8„) output of LANDSAT. 
Notice the heavy, regular striping. 

Destriped image of band 6 output. 

The destriping lookup tables for band 6 - inverse transfer functions 
for the six image sensors. The nonlinearities of the transducers are 


apparent. 


Original, unrectified image of band 7 (.Bp to ].1„) output of LANDSAT 
Notice regular striping. 

Destriped image of band 7 output. 

The destriping lookup tables for band 7 -- inverse transfer functions 
for the six image sensors. These transducers appear to be fairly 


linear, differing mostly only as 


regards gain or amplification. 
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FIGURE 2 









































